#!/usr/bin/perl
# 2014/06/20 To display association in GenomeBrowser, convert map to BED for specific region

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 4 or $ARGV[0] eq '--help'){
	die "Usage : $0 -t [target specified file] -i [map file]\n";
}

my %opt;
getopts("i:t:", \%opt);
my $FILE_target = $opt{t};
my $FILE_map = $opt{i};



#---------------------------------------
# read target file
#---------------------------------------
my @target;
my $fh_target = IO::File->new($FILE_target) or die "cannot open $FILE_target: $!";
while($_ = $fh_target->getline()){
	s/\r?\n//;
	my ($name, $chr1, $start1, $end1, $chr2, $start2, $end2) = split /\t/;
	if($chr1 eq $chr2){
		if($start1 < $start2){
			push @target, [$chr1, $start1, $end1, $chr2, $start2, $end2];
		}else{
			push @target, [$chr2, $start2, $end2, $chr1, $start1, $end1];
		}
	}else{
		if($chr1 eq "I" or ($chr1 eq "II" and $chr2 eq "III")){
			push @target, [$chr1, $start1, $end1, $chr2, $start2, $end2];
		}else{
			push @target, [$chr2, $start2, $end2, $chr1, $start1, $end1];
		}
	}
}
$fh_target->close();


#---------------------------------------
# extract target and output as BED format
#---------------------------------------

my $fh_map = IO::File->new($FILE_map) or die "cannot open $FILE_map: $!";
while($_ = $fh_map->getline()){
	s/\r?\n//;
	my ($id, $chr1, $loc1, $direction1, $mapQ1, $hinID1, $hinLoc1, $chr2, $loc2, $direction2, $mapQ2, $hinID2, $hinLoc2, $type, $num, $locations) = split /\t/;

	foreach my $t(@target){
		my ($t_chr1, $t_start1, $t_end1, $t_chr2, $t_start2, $t_end2) = @$t;
		if($t_chr1 eq $chr1 and $t_start1 < $loc1 and $loc1 < $t_end1 and $t_chr2 eq $chr2 and $t_start2 < $loc2 and $loc2 < $t_end2){
			my $second_block_start = $loc2 - 100;
			print join("\t", $chr1, $loc1, $loc2, $id, 0, "+", $loc1, $loc2, "255,0,0", 2, "100,100",  "0,$second_block_start") . "\n";
		}elsif($t_chr1 eq $chr2 and $t_start1 < $loc2 and $loc2 < $t_end1 and $t_chr2 eq $chr1 and $t_start2 < $loc1 and $loc1 < $t_end2){
			my $second_block_start = $loc1 - 100;
			print join("\t", $chr1, $loc2, $loc1, $id, 0, "+", $loc2, $loc1, "255,0,0", 2, "100,100",  "0,$second_block_start") . "\n";
		}

#		if($t_chr1 eq $chr1 and $t_start1 - 5000 < $loc1 and $loc1 < $t_end2 + 5000 and $t_chr2 eq $chr2 and $t_start1 - 5000 < $loc2 and $loc2 < $t_end2 + 5000){
#			my $second_block_start = $loc2 - 100;
#			print join("\t", $chr1, $loc1, $loc2, $id, 0, "+", $loc1, $loc2, "255,0,0", 2, "100,100",  "0,$second_block_start") . "\n";
#		}elsif($t_chr1 eq $chr2 and $t_start2 - 5000 < $loc2 and $loc2 < $t_end1 + 5000 and $t_chr2 eq $chr1 and $t_start2 - 5000 < $loc1 and $loc1 < $t_end1 + 5000){
#			my $second_block_start = $loc1 - 100;
#			print join("\t", $chr1, $loc2, $loc1, $id, 0, "+", $loc2, $loc1, "255,0,0", 2, "100,100",  "0,$second_block_start") . "\n";
#		}


	}
}
$fh_map->close();